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SUMMARY 

The paper presents a numerical method useful to describe unsteady three- 
dimensional flow fields within turbomachinery stages. The method solves the 
compressible, time-dependent, Euler conservation equations with a finite-volume, 
f lux-splitting, total-variation-diminishing, approximately factored, implicit scheme 
Multiblock composite gridding is used to partition the flow field into a specified 
arrangement of blocks with static and dynamic interfaces. The code is optimized to 
take full advantage of the processing power and speed of the Cray Y/MP supercomputer 
The method is applied to the computation of the flow field within a single-stage, 
axial flow fan, thus reproducing the unsteady three-dimensional rotor-stator 
interaction. 


NOMENCLATURE 

A,5,C quasi-linear matrices 
a^ strength of wave 

b^^^ speed of wave 

E total internal energy 

e internal energy 

F,G,H flux vectors 
f unknown vector 

H total enthalpy 

h enthalpy 

IB leading edge i index 

IE trailing edge i index 

i,j,k space discretization indices 
J Jacobian 

L^ limiter 

left eigenvector 
M Mach number 

m plane index 
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NI maximum block i count 

NJ maximum block j count 

NK maximum block k count 

p pressure 

right eigenvector 
T residual vector 

t time 

U,V,W contravariant velocity components 
UfVfW velocity components 
x,y ,2 Cartesian coordinates 

a specific heat ratio 

limiter parameters 
r^n^n curvilinear coordinates 

5t time step 

5x mesh size 

8,^ time discretization parameter 

summation over positive wave speeds 
E" summation over negative wave speeds 

(7 density 

$ space discretization parameter 

Subscripts ; 

1 left 

r right 

Superscript : 
n time level 


INTRODUCTION 

The numerical methods for the simulation of single-domain, steady, three- 
dimensional flow fields developed over the past years are now widely used within a 
design environment. These methods include both Euler (Denton 1983; Van Hove 1983) 
and Navier-Stokes (Dawes 1986 a,b, 1987; Hah 1984, 1986; Shamroth et al. 1982; Roscoe 
et al . 1984; Weinberg et al. 1986; and Moore 1985) solvers. Although the Navier- 
Stokes methods can produce better computational results, their use as a design tool 
appears to be certainly more difficult, due to the increased number of parameters 
influencing the solution accuracy. Analyses of flows through isolated rows can be 
used to study many of the aerodynamic phenomena in turbomachinery, but such analyses 
yield no information regarding the unsteadiness arising from the interaction of 
rotating and stationary blade rows. 


2 


Since the flow within a turbomachine is extremely complex, computational methods 
can be used to great advantage in understanding such flows. An accurate analysis of 
the flow associated with rotating and stationary row configurations can be very help- 
ful in optimizing the performance of the whole turbomachine. Within each stage, 
rotor and stator rows alternate and thus produce particularly strong interactions. 

As shown by Dring et al. (1982), the temporal pressure fluctuation near the leading 
edge of the rotor can be as much as 72 percent of the exit dynamic pressure in a 
rotor-stator interaction problem with a 15-percent chord length axial gap. Numerical 
methods to simulate multidomain, unsteady, three-dimensional flow fields appear 
therefore to be particularly interesting. The present paper describes the develop- 
ment of one of these new numerical methods. 

The simulation of the rotor-stator interaction problem needs the introduction of 
a particularly complex flow model. The problem is assumed to be governed by the 
three-dimensional Euler conservation equations, written in unsteady, compressible 
form, with the gas density, velocity, and energy as the basic variables. These equa- 
tions have to be solved in a particularly complex flow domain. Multi-blade-passage 
and multi-blade-row configurations have to be considered to simulate the unsteady 
rotor-stator interaction. Turbomachinery bladings are usually highly skewed, and the 
hub and tip surfaces strongly diverge. Moreover, fixed and movable blade passages 
alternate. As a result, both fixed and movable blade passages are bounded by 
strongly irregular surfaces, and information must be transferred from a fixed to a 
movable reference frame. Furthermore, high flow resolution in both space and time is 
required. Finally, the unsteady flow behavior requires the proper treatment of 
boundary conditions. 

The development of a numerical method resolving the flow model previously 
described poses many programming problems. The equations have to be solved in a 
time-dependent, body-fitted, curvilinear, reference frame. The discretization of a 
flow domain including multi-blade-passages and multi-blade-rows requires the intro- 
duction of multiblock gridding, where the whole flow domain is partitioned in a 
specified arrangement of blocks. These blocks are limited by the boundaries of the 
flow domain or by block interfaces. The high resolution in space requires, in addi- 
tion to the introduction of a proper grid refinement where higher flow gradients are 
expected to occur, a proper selection of block arrangements and individual block 
dimensions. The use of a finite-volume space discretization, with an accurate flux- 
corrected interface flux-splitting formula and a total-variation-diminishing approach 
in limiting the components of the interface flux, produces high spatial accuracy. A 
multistep, approximately factored, implicit scheme produces high time accuracy. 
Multiblock gridding with static and dynamic grid interfaces is used to properly 
discretize the complex flow field. The relative motion of rotor and stator blades is 
simulated by using grids that move relative to each other, with proper treatment at 
the interface between blade rows. The implementation of appropriate boundary condi- 
tions is based on characteristic concepts. A characteristic variable boundary condi- 
tion technique is introduced in the time-dependent, body-fitted reference frame for 
inflow, outflow, and solid boundaries, and it is implemented by using phantom cells. 

The method is used to model the unsteady flow in a single-stage, transonic, 
axial flow fan. The test case considered here does not exhibit such a particularly 
strong interaction as that shown in the work of Dring et al . (1982). The rotor and 
stator are separated by approximately 85-percent rotor chord at midspan, thus reduc- 
ing (but not eliminating) the flow field interactions between blade rows. The compu- 
tational domain is made up of two rotor blade passages, each discretized by using 
about 13 500 grid points, and three stator blade passages, each discretized by using 
8500 grid points. The total number of grid points involved in the rotor-stator 
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problem discretization therefore equals 52 500. Approximately 190 vital pieces of 
information have to be known for each cell. The computations have been performed on 
the NASA Ames Cray Y/MP computer. At the beginning of the computations, the rotor- 
stator interface orientation is defined with the rotor-blade trailing edges located 
midway between the neighboring stator— blade leading edges. The user CPU time 
required for each rotor-stator orientation is about 700 sec. The system CPU time is 
about 150 sec. The internal memory required is only 2 megawords. The result 
obtained after four interface reorientations are presented and compared with laser 
anemometer measurements for both rotor (Pierzga and Wood 1985) and stator (Suder 
et al, 1987/ and Hathaway et al. 1987) rows. The theoretical results are also 
compared with other steady-state theoretical results. These comparisons allow a 
first assessment of the code prediction capability. 


EULER FLOW MODEL 

The computation of high Reynolds number, unseparated, unsteady cascade flow 
fields can easily be accomplished by solving the Euler conservation equations, 
written in unsteady, compressible form. The basic variables are the gas density, 
velocity, and energy. Their conservation equations are written in a body-fitted, 
time-dependent, curvilinear, reference frame: 


r = r (x,y, z, t ) 
fl «* [l{x,y, z, t ) 

n = ri(x,y,z,t) 


In vector form, 
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h = H - (u^ + v^ t w^)/2 
e = E - (u^ + v^ + w^)/2 

The Jacobian J, the metric quantities T T F 0^^, £1 0^ 

n , and n , and the contravariant velocities U, V, and W are given in 

appendix A. 


The integration requires the introduction of the conditions to be fulfilled 
along the boundaries of the domain of interest. The physical flow domain is limited 
by inflow, outflow, solid, and periodic boundaries where appropriate boundary condi- 
tions are needed. The boundary conditions along inflow, outflow, and solid bound- 
aries are defined according to characteristic concepts - i.e., by discretizing the 
flow equations written in characteristic form, thus simulating the correct unsteady 
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flow behavior. Before starting the time integration, proper initial conditions have 
to be specified. These conditions can be either a previously computed steady-state 
condition or a rough approximation of a steady flow condition. 


NUMERICAL METHOD 

The integration of the Euler equations, over discrete contiguous volumes, in a 
computational space where 5 x=5y=5z= 1, yields for the grid point i,j,k the follow- 
ing cell-centered finite-volume formula; 


f = 

-5 P - 

5 G - ( 

/t 

where 

1 

j 

5 F = 

F 

- F 

i 

i+l/2 

i-l/2 

5 G = 

G 

- G 

j 

j + l/2 

j-l/2 

5^h = 

®k+l/2 

- »K-l/2 


The dependent variable vector f represents values at the cell centers, while the 
vector flux functions F(f), G(f), and H(f) represent values at the interface 
between neighboring cells. The evaluation of the vector flux functions at the cell 
interfaces is performed by using the values in the neighboring cells. The interface 
terms are computed according to a one-parameter flux-difference-splitting scheme. 

For ease of understanding, the flux-difference-splitting scheme is presented 
here in one space dimension. As reported in appendix B, a family of high spatial 
accuracy formulae is written as follows: 


pi A ss F 

i+1/2 i+1/2 


+ (1 + + (1 - ♦)/4(dF\.^/2 - dF-^^3^^ 


where dF is the interface flux difference for the collection of waves. The super- 
script denotes either positive or negative traveling waves. The option $=l/3 gives 
a third-order accurate scheme; the option $=-l gives a fully upwind second-order 
accurate scheme. These one-dimensional formulae can be extended to three dimensions 
by assuming that all waves travel normal to their respective interfaces. 

Spurious oscillations are controlled by using total-variation-diminishing con- 
cepts in limiting the components of the interface flux. The use of limiters yields 
the following expressions for the corrective flux terms; 

dF' 

dF' 

dF' 

dF' 

It J/ z 

where is a right eigenvector of the quasi-linear matrix A=F Different 

limiter expressions have been proposed. These expressions are presented in 
appendix C. 
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The previous equations are finally integrated in time by using an approximately 
factored, implicit scheme. The equation is first written in the linearized, dis- 
crete, integral delta form 

(I + 05t/(l + (/>)»"*]■ 5f" = -5t/(l + ^)T" + f/(l + 

Different members of this family have particular relevance* If 0=1 and ^=1/2, the 
scheme is three-point backward; if 0=1 and ^=0, it is backward Euler; and if 0=1/2 
and (p=0f it is trapezoidal. The left side of the previous time-stepping formula is 
expressed as follows: 


m""* = d n** t 5 A"* + 6b** + 5 B"* + + 5 c’* 

i 1 1 j 


where 
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with \ and resulting from flux-vector-splitting theory {Steger and 

Warming, 1981) while the right side is expressed as 


t"" = S.f’" + 5,g" + 
i j 

with F, G, and H resulting from the flux-difference-splitting theory discussed 
herein (for further details, see Janus 1989). 


For ease of computer programming, the left side operator is split into the 
product of two operators: 

(I + 05t/<i + = -5t/(i + ^)t" + ^/{i + 

(I + 05t/(i + ^)M'*]-6f" = 5f" 

= f" + 5f" 


In the solution procedure, since the first operator is applied by moving forward, it 
is referred to as the forward operator. The second one is applied by moving backward 
and therefore is referred to as the backward operator. The left side operators are 
defined as follows: 

+ (5 c*'’’* 

1 j k 

Despite factoring, the scheme seems to retain the original unconditional stability. 

The implementation of appropriate boundary conditions is based on characteristic 
concepts. A characteristic variable boundary condition technique is introduced in a 
time-dependent body-fitted reference frame for inflow, outflow, and solid boundaries 
(Whitfield and Janus 1989). These boundary conditions are implemented by using 
phantom cells, with changes in dependent variables and 5f* to be set equal to 

zero in these phantom cells. It is pointed out that the outflow boundary condition 
uses both the characteristic variable form of the equations and a radial equilibrium 
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condition where the outflow static pressure at the hub is specified. Better details 
on the numerical method can be found in Janus (1989). 


COMPUTER CODE 

The code was written (Whitfield and Janus 1989) and further optimized for use on 
the NASA Ames Cray Y/MP computer. The Cray Y/MP has 8 processors, 32 megawords of 
common memory and 256 megawords of solid-state storage, and it allows particularly 
efficient computational fluid-dynamic simulations, especially when some guidelines 
are followed in the code development. 

Because of the in-core memory limitations associated with the computation of 
complex three-dimensional flow fields, the overall flow environment is segmented into 
several smaller and more manageable intercommunicating flow environments. Figure 1 
shows the composite field circumferential and axial partitioning. The global domain 
is first circumferentially partitioned by considering different blade passages. The 
domain is then axial partitioned by considering different blade rows. Each blade row 
is granted limited rotational freedom relative to the adjacent blade rows. When the 
domain is completely partitioned, a block referencing scheme is introduced following 
the global axes orientation. Within each block, the local axes are assumed to follow 
the global axes. The block index limits NI, NJ, and NK are restricted such that 
for all blocks within a blade row NI, NJ, and NK remain constant. In addition, NJ 
remains constant between blade rows. Variations in NI are permitted between blade 
rows. Variations in NK are permitted between blade rows in the limits of a con- 
stant circumferential cell count (the product of NK by the number of circumferen- 
tial blocks) to be satisfied. 

The reduction of in-core (primary) memory requirements calls for the use of 
secondary memory, in particular for the use of the Cray's rapid access solid-state 
storage device (SSD). The use of computational blocks of differing cell count, and 
therefore differing memory requirements, complicates memory utilization. The ribbon 
vector dynamic memory management (Janus 1989) stores the entire field on secondary 
memory with the exception of the block currently under execution. The method allows 
the adjustment of the in-core field length to accommodate each block, regardless of 
size, without excess memory words. 

The optimization of the processing of the data residing in memory calls for a 
reduction of 10 wait time. Unblocked data transmission to and from secondary memory 
is used. All data contained within the ribbon vector are transferred between in-core 
and secondary memory via unblocked standard FORTRAN I/O statements. Proper unblocked 
data length has to be provided. 

Vectorizing FORTRAN compilers identify inner DO loops that are suitable for 
vectorizat ion . Since this enables the loop to run in vector mode, the program runs 
much faster. It is important to note that if a loop has even one nonvector izble 
construct the entire loop will not be vectorized. There are a number of constructs 
that can decrease the performance of a code. The process of optimization of a code 
to obtain peak performance on the Cray Y/MP can be regarded as the elimination of 
these constructs. Some general vectorizing guidelines are provided for this purpose. 
The constructs that can decrease the performance of a code are usually referred to as 
memory strides, vector lengths, vector dependencies and recursions, logically sepa- 
rate vectors within a single array, loops with IF statements, loops with subroutines, 
and functions calls . 
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The term stride refers to the increment in memory between successive words 
fetched or stored in a vector operation. Bad strides can cause the code to run sig- 
nificantly slower. The vector length of an operation is the number of times an inner 
DO loop is executed. Longer vectors are processed more efficiently. DO loops cannot 
be vectorized in cases where calculations for one iteration of a loop require results 
from a previous iteration (vector dependency). Loops that reference independent 
sections of a single array can often appear to have vector dependencies. IF state- 
ments prevent loops from being vectorized* Subroutine and function calls prevent 
vectorization for many reasons. 

The code development has to include some general guidelines. Vectors have to be 
thought of as fundamental constructs. Vector operations have to be isolated from 
scalar operations. Control statements have to be avoided in loops. Recursion has to 
be avoided or moved to outer loops. Power-of-two memory strides have to be avoided. 
The array subscripts have to be kept simple and explicit. In nested DO loops, the 
inner DO loop has to be made as long as possible. Double precision has to be avoided 
unless absolutely necessary. Floating-point arithmetic operations have to be used 
instead of integer arithmetic operations (except in array subscripts). System 
library routines have to be used when possible. 

In the cases of loops that reference independent sections of a single array, 
indexing techniques should be used to define two logical arrays within a single 
physical array. There are many approaches to vectorizing loops with IF statements. 
The method called loop splitting divides a loop into its vectorizable and nonvector- 
izable parts. The loop restructuring removes redundant IF statements or restructures 
loops with invariant IF statement. The approaches to vectorizing loops with sub- 
routine calls include using statement functions, promoting the subroutine into the 
loop, or expanding the subroutine to include the loop. 

The adopted approximately factored, implicit scheme poses a problem of vector 
dependency. The implicit scheme involves point simultaneous solutions and backward 
or forward substitutions, and it poses a problem of vector dependency during the 
substitution. To circumvent this, a procedure simultaneously processing those cells 
lying on a special diagonal plane in the computational space is introduced. Cells 
whose indices satisfy the equation of the diagonal plane, expressed as i+j+k=m, where 
m is a constant designating the plane level, are computationally independent. 
Diagonal plane processing can then be used with the aid of indirect addressing to 
facilitate the vectorization of a backward or forward substitution. 

The advance of the flow field in time requires two global sweeps through the 
global domain, one forward and one backward. In addition, all blocks are forward 
swept internally and then they are backward swept internally. An equivalence to the 
analogous single-block solution is thus provided. The internal forward sweep within 
each block consists of operating on each cell on a diagonal plane, from the lowest 
level interior plane to the highest level interior plane, by applying the forward 
operator at each computational cell. After information due to communicat ion has been 
collected, the in-core resident block is returned to secondary memory and the next 
block is transferred in-core. The internal backward sweep is similar, although the 
planes are traversing in decreasing rather than in increasing order, and the operator 
applied is the backward one. 

The communication between blocks is accomplished by considering continuous grid 
lines between blocks. Values from within the domain of one block are extracted and 
then injected as phantom values in an adjacent block. Since each blade row can 
rotate relative to its adjacent blade rows, continuous grid lines are required across 
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the shearing block-block interface at all time levels. The shearing block-block 
interface in physical space yields a dynamic interface in computational space, with a 
progressive changing of communication partners. The block-block communication is 
properly accomplished by simulating the internal cell communication within the 
blocks. The implementation of boundary conditions as well as the transfer between 
neighboring blocks is implemented by using phantom cells. This approach wastes a 
small amount of memory for the offset, fully justified by the simplification of the 
problem. Other details about the computer codes can be found in Janus (1989). 


RESULTS 

The code has been applied to the computation of the flow in a single-stage, 
transonic, axial-flow fan. Figure 2 shows the plane view of the single-stage fan 
flow path (Hathaway 1986). Figure 3 shows the representative rotor blade sections at 
three spanwise locations, hub, midspan, and tip. Figure 4 shows the representative 
stator blade sections at three spanwise locations, hub, midspan, and tip. The design 
speed of the fan is 16 043 rpm and the corrected mass flow rate is 34 kg/s. The 
rotor and stator are separated by approximately 85-percent rotor chord at midspan. 

The advantage of the wide axial spacing between blade rows is that the flow field 
interactions between the blade rows are reduced, thus decreasing the fan noise. 

The rotor is composed of 22 blades of multiple-circular-arc design. The rotor 
aspect ratio (averaged) is 1.550. The inlet tip diameter is 51.3 cm; the inlet hub/ 
tip radius ratio is 0.375; and the rotor tip clearance is 0.5 mm. The stator is 
composed of 34 blades of double-circular-arc design. The stator discharges the fluid 
axially. The axial chord is a constant 5.6 cm from hub to tip. The tip diameter is 
constant at 48.7 cm. The inlet hub/tip radius ratio is 0.500, while the exit hub/tip 
radius ratio is 0.530. 

Some unsteady experimental results (Hathaway et al. 1987) are presented in 
figure 5, The dominant blade-row interactions are due to viscous interactions caused 
by the chopping of rotor wakes by the downstream stator blade row. The figure shows 
a sequence of plots of the turbulence kinetic energy contours, for a rotor rotation 
of one pitch, at midspan, SP denotes rotor shaft position: 50 rotor shaft positions 

are specified per rotor pitch. The shaded regions identify rotor wake fluid. After 
the stator blade chops the rotor wake, the wake segments move at different speeds 
along the stator blade pressure and suction sides. When the rotor wake segments 
reach the stator exit, there is a drift between the rotor wake segments originally 
part of the same rotor wake. 

Calculations have been performed by considering two rotor and three stator blade 
row blocks. Each rotor blade row block is made up of NixNJxNK grid points, with 
NI=49 planes from inlet to exit, NJ=21 planes from hub to tip, and NK=13 planes from 
suction to pressure side. The leading edge is located at IB=21, while the trailing 
edge is located at IE=41, Each stator blade row block is made up of NI=45 planes 
from inlet to exit, NJ=21 planes from hub to tip, and NK=9 planes from suction to 
pressure side. The leading edge is located at IB=9, while the trailing edge is 
located at IE=21, The last five rotor axial planes and the first five stator axial 
planes are used for communication between blade rows. The rotor and stator blade 
passage computational grids are shown in figure 6. 

Calculations are started from a rough approximation of the flow field and the 
flow equations are then integrated in time. The selected code options allow second- 
order accuracy in time and third-order accuracy in space. The limiter adopted is the 
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van Leer limiter* The relative motion between blade rows requires 264 time cycles 
for each interface reorientation- The computations have been performed on the NASA 
Ames Cray Y/MP computer. The user CPU time required for each interface reorientation 
is about 700 sec, while the system CPU time is about 150 sec. The internal memory 
required in only 2 megawords. 

The upstream boundary conditions are specified in terms of a free stream Mach 
number, taken equal to 0.5565, and free stream flow angles, defining an inlet axial 
direction. The downstream boundary condition is specified in terms of a hub static 
pressure, fixed at 1.7869 times the inlet static pressure. These boundary conditions 
simulate the real flow conditions with some difficulties. For many reasons, the 
correct mass flow rate is not perfectly verified. The mass flow rates experimentally 
measured have a relatively strong degree of uncertainty. Many runs are required to 
establish the proper boundary conditions to obtain a prefixed value in the mass flow. 
Most of the results available for the rotor have been obtained without a downstream 
stator row. Finally, it is well known that the best agreement with experimental data 
in computing flows within compressors or fans requires the evaluation of viscous 
effects (Pierzga and Wood 1985). For all these reasons, the boundary conditions have 
been defined on the basis of simple one-dimensional concepts, without running the 
code, and the correct mass flow rate is therefore not verified in the calculations. 

The three-dimensional rotor-stator interaction simulation provides the user a 
huge amount of information. In a first stage of validation, only the results 
obtained for a particular rotor shaft position are considered. The results obtained 
after four interface reorientations are presented and compared with laser anemometer 
measurements for both rotor (Pierzga and Wood 1985) and stator rows (Hathaway et al. 
1987). These comparisons allow the user to assess the prediction capability of the 
Euler code not only from a qualitative but from a quantitative viewpoint. 

Figure 7 to 16 show some computational results obtained for the rotor row, com- 
pared with available experimental data and other computational results, obtained 
using Denton's code (Denton 1983). The experimental results were obtained with a 
laser anemometer for a rotor row without a downstream stator row (Pierzga and Wood 
1985). Denton's Euler code, modified to include the effects of boundary layer dis- 
placement (in order to obtain the best agreement with experimental data), was applied 
to the test case defined in the experiments (Pierzga and Wood 1985). 

Figure 7 shows (from left to right) the relative Mach number contours in the 
rotor row leading edge, trailing edge, and exit planes. The rotor flow appears to be 
strongly three-dimensional and rotational. 

Figure 8 shows the relative Mach number contour results obtained with the laser 
anemometer (left) and Denton's code (right) at 30-percent span from the tip. Fig- 
ure 9 shows the present computational results for the same flow surface. The rela- 
tive inlet Mach number is supersonic. The shock location is shown to provide a 
clearer picture of the flow in the passage. The shock location was determined by 
considering the Mach number and flow angle data in the streamwise direction and by 
assuming the starting point of the flow deceleration as the shock front (Pierzga and 
Wood 1985). The agreement between data is quite good, with only minor differences. 

A normal shock is followed by a second shock. The peak Mach number ahead of the 
first shock is about 1.4. This shock is accurately located by the present code with 
a peak Mach number of about 1.375. The second shock is shown extending from the 
pressure surface to about midpitch, while closer to the suction surface, the location 
of the shock is more difficult to determine. The location and the extension of this 
second shock is the major difference between the experiment and the analysis. The 
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computed Mach number in the rear part of the blade is only slightly lower than the 
measured one, 0.85 compared to 0.90. The numerical calculations produce shocks that 
are smeared over several grid points, probably due to grid coarseness and neglected 
viscous effects. The code prediction capability appears to be quite good, especially 
when compared with the accuracy provided by Denton's code, certainly more calibrated 
on the specific test case. There is a higher static pressure downstream of the rotor 
blade row that reduces the flow velocities in the rear part of the blade passage. 

A more quantitative comparison between experimental and theoretical data is the 
blade surface velocity comparison. The blade surface velocity data predictions are 
the most important in helping the designer tailor the blade shape. Figure 10 shows 
the comparison between experimental and theoretical results for the relative Mach 
number differences between pressure and suction side at 30-percent span from the tip. 
The laser anemometer data are taken at the first point off the blade surface, close 
to 5- and 95-percent pitch, and therefore the computational results are selected at 
the same locations. The agreement appears to be quite satisfactory in both the 
suction and pressure side velocities, with only minor changes in the values on the 
rear blade surfaces. 

Figure 11 shows the relative Mach number contour results obtained with the laser 
anemometer (left) and Denton's code (right) at 10-percent span from the tip. Fig- 
ure 12 shows the present computational results for the same flow surface. The rela- 
tive inlet Mach number is increasingly supersonic with reference to the value in the 
section at 30-percent span from the tip. The agreement between data is particularly 
good. The presence of a second shock is not clearly evident. The peak Mach number 
ahead of the first shock is about 1.4. The shock is accurately located by the 
present code with a peak Mach number of about 1.375. The computed Mach number in the 
rear part of the blade is now lower than the measured one, 0.875 compared to 0.975, 
but the code prediction capability appears to be quite satisfactory, within the 
limits of accuracy obtained in a solution of the Euler equations. 

Figure 13 shows a comparison between experimental and theoretical results for 
the relative Mach number differences between pressure and suction side at 10-percent 
span from the tip. The agreement appears to be better than that obtained on the flow 
surface at 30-percent span from the tip. Minor differences in the values on the rear 
blade surfaces still remain. 

Figure 14 shows the relative Mach number contour results obtained with the laser 
anemometer (left) and Denton's code (right) at 70-percent span from the tip. Figure 
15 shows the present computational results for the same flow surface. The relative 
inlet Mach number is now subsonic. The peak Mach number close to the blade leading 
edge is lower than that in the experiments, as well as the rear flow Mach number, 
especially close to the blade pressure side. The mass flow close to the hub seems to 
be strongly reduced with reference to the value obtained in the experiments. 

Finally, figure 16 shows the comparison between experimental and theoretical results 
for the relative Mach number differences between pressure and suction sides at 
70-percent span from the tip. The agreement appears to be worse than that obtained 
on the flow surfaces at 10- and 30-percent spans from the tip, and there are strong 
differences in surface velocities especially on the blade pressure side. 

Most of the differences between experimental and computational results seem to 
be due to the influence of the downstream stator row on the rotor flow field. The 
flow picture provided by the Euler flow model certainly appears to be better than 
expected. 
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While experimental data are available in enough detail for the rotor flow (but 
these results do not consider a downstream stator row), only a few experimental data 
have been provided for the stator row. Despite this, these results have a particular 
significance since they have been obtained by considering an upstream rotor row. 

Figures 17 to 19 show some computational results obtained for the stator row. 
Figure 17 shows (from left to right) the Mach number contours in the stator row 
inlet, leading edge, and trailing edge planes. The stator flow appears to be 
strongly three-dimensional and rotational from the inlet plane. Figure 18 shows the 
Mach number contour results obtained with the laser anemometer (Hathaway et al, 1987) 
at 50— percent span from the tip. Figure 19 shows the present computational results 
for the same flow surface. The inlet Mach number is subsonic. The agreement between 
data is quite good. The peak Mach number on the suction side is predicted fairly 
well, with a very similar 0.70 constant Mach number contour. The same is true for 
the Mach number in the rear part of the blade, with only minor underestimation, 0.55 
compared to 0.59. 

The relative influence of stator and rotor rows appears to be well simulated. 


CONCLUSIONS 

This paper has presented a computer code developed for the analysis of unsteady 
three-dimensional flow fields within turbomachine stages. The model takes into 
account the unsteady flow fields within complex domains including multi-blade-passage 
and multi-blade-row (both fixed and movable) configurations. The flow model has a 
particularly wide generality and is thus applicable to the solution of the majority 
of problems arising in the aerodynamic and acoustic design of turbomachinery 
components . 

The multiblock gridding allows one to properly discretize the multi-blade- 
passages and multi-blade-row configurations. The partitioning of the whole flow 
domain in a specified arrangement of blocks is the only answer to discretize flow 
domains requiring a huge number of grid points with strongly varying grid refine- 
ments. The block interface treatment probably requires further work to improve the 
transfer of information without introducing interface errors. 

The finite-volume, flux-corrected interface flux-splitting, total-variation- 
diminishing space discretization allows up to third-order accurate space discretiza- 
tions without undesired artificial viscosity. The multistep, approximately factored, 
implicit time discretization shows an apparently unconditional stability and produces 
up to second-order accuracy in time. Finally, the efficient use of the hardware 
capability of the Cray y/MP supercomputer leads to significantly reduced costs of 
very complex fluid-dynamic simulation. 
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APPENDIX A - COORDINATE TRANSFORMATION 


The Euler conservation equations are written in a body-fitted, time-dependent, 
curvilinear, reference frame. The coordinate transformation 

r = r (x,y, z,t) 
t) - Q(x,y,z,t) 
n = ri(x,y, z,t) 

leads to the following Jacobian and metric quantities expression; 
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The contravariant velocities are given as follows: 
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APPENDIX B - PLUX-DIFFERENCE-SPLITTING FORMULAE 


The evaluation of the vector flux functions at the cell interfaces is performed 
by using the values in the neighboring cells according to a one-parameter flux- 
difference-splitting scheme. For ease of understanding, the scheme is presented here 
in one space dimension. The flux difference at the interface between left and right 
cells is expressed as 


dF = F - F = A{f - f,) = A-df 

where A=F ^ is s quasi-linear matrix, representative of local interface. The 
matrix A ' is evaluated at the interface according to the Roe averaging procedure 
(Roe et al. 1984) by using the following relations; 


a 

u 

H 


+ ay\)/(a 


1/2 

1 

1/2 

1 


+ 

+ a 


1/2 

r 

1/2 


) 

) 


Each set of left or right eigenvectors forms a spanning set in state space, and a 
basis is constructed using these orthonormal vectors. The interface differential df 
is thus proportional to the right eigenvectors r of A; 

df = E 

where a^ is the magnitude of the component in the direction (strength of the 

wave, i.e., the jump in the characteristic variable across the interface). The 
interface flux difference can therefore be expressed as the composition of a collec- 
tion of waves as follows: 

dF = Ea^b^^^r^^^ = dF^ + dF' = E^a^b^^^r^^^ + E”a^b^^ ( j ) 

where is an eigenvalue of A (speed of j^^ wave). Symbols E^ and E denote 

summation over positive and negative wave speeds, respectively. The interface flux 
can therefore be computed from one of the following first-order expressions: 

F = F + E”a 

i+l/2 1 j 

F ,,, = F - E*a 

i+l/2 r j 

= 1/2 [f + F - Ea lb‘^’ lr‘^’1 

i + l/2 'Ll r j ' ' 

The addition of a corrective flux to the previous formulae produces higher 
spatial accuracy flux-difference-splitting formulae. The family of these formulae is 
written as follows; 


= ^.1/2 + <1 ^ - dF-^,,/3) ^ ^ *)/ 4 (dFVl /2 

where the principal part of the truncation error is 
(1/3 - ^)/4{5x)^r ^ 


dF 


i-+3/2 
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APPENDIX C - FLUX LIMITERS 


Total-variation-diminishing flux expressions make use of flux limiters* These 
limiters are generally made dependent on the values of a parameter e proportional 
to the change in the characteristic variables across nearby interfaces: 


.(J) 


IT'-' = a = b 

i +*/2 l +*/2 J l +*/2 


= 1 

i +*/2 




*df 

i +*2 i + y */2 


where ] 
follows J 


iU) 


is the left eigenvector of A, The "minmod” 


limiter is defined as 


minmod (x,y) = sign(x)max{0, min[|x|, ysign(x)]} 

where P is a compression parameter assumed to be equal to (3 - $)/(! - $)• The 
"superbee" limiter is defined as follows: 

L^(tn,n) = 

cmplim(x,y) = sign ( x) max{0 , min[|x|, ^ysign(x)], min[^|x|, ysign(x)]} 

where the compression parameter is now assumed to be ^=2. The "van Leer" limiter is 
finally defined as follows: 

L^(tn,n) = 

vanlim(x^y) = (xy + |xy|)/(x + y) 
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Figure 1. — Axial and circumferential partition of the flow domain. 
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Rgure 3. — Rotor blade sections. 



Figure 4. — Stator blade sections. 
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Figure 9. — Theoretical (present code) relative Mach number contours in rotor 
30-percent span from tip (results obtained with downstream stator row). 
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Figure 10. — Comparison between theoretical (present 
code, results obtained with a downstream stator row) 
and experimental (Pierzga and Wood, 1985, results 
obtained without a downstream stator row) blade 
relative Mach number distributions in the rotor row at 
30-percent span from the trip. 
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Figure 14.— Experimental (left) and theoretical (Denton’s code, right) relative Mach number contours in rotor row at 70-percent span from 
tip (Pierzga and Wood, 1985, results obtained without downstream stator row). 



Figure 15. — Theoretical (present code) relative Mach number con- 
tours in rotor row at 70-percent span from tip (results obtained 
with downstream stator row). 



Figure 16. — Comparison between theoretical (present 
code, results obtained with a downstream stator row) 
and experimental (Pierzga and Wood, 1985, results 
obtain^ without a downstream stator row) blade 
relative Mach number distributions in the rotor row at 
70-percent span from the trip. 





Figure 18. — Experimental Mach number coniours in stator row at 
50‘percent span from tip (Hathaway et al., 1987, results obtained 
with an upstream rotor row). 
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